A framework for large-scale relativistic simulations in the characteristic approach 
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We present a new computational framework (LEO), that enables us to carry out the very first 
large-scale, high-resolution computations in the context of the characteristic approach in numerical 
relativity. At the analytic level, our approach is based on a new implementation of the "eth" 
formalism, using a non-standard representation of the spin-raising and lowering angular operators 
in terms of non-conformal coordinates on the sphere; we couple this formalism to a partially first- 
order reduction (in the angular variables) of the Einstein equations. The numerical implementation 
of our approach supplies the basic building blocks for a highly parallel, easily extensible numerical 
code. We demonstrate the adaptability and excellent scaling of our numerical code by solving, 
within our numerical framework, for a scalar field minimally coupled to gravity (the Einstein-Klein- 
Gordon problem) in the 3-dimensions. The nonlinear code is globally second-order convergent, and 
has been extensively tested using as reference a calibrated code with the same boundary-initial data 
and radial marching algorithm. In this context, we show how accurately we can follow quasi-normal 
mode ringing. In the linear regime, we show energy conservation for a number of initial data sets 
with varying angular structure. A striking result that arises in this context is the saturation of the 
flow of energy through the Schwarzschild radius. As a final calibration check we perform a large 
simulation with resolution never achieved before. 

PACS numbers: 04.25.Dm, 04.30.Db, 04.70.Bw 04.20.Ex, 04.25. Nx 

Keywords: black holes, characteristic evolution, spherical coordinates, finite difference methods, parallel 

computation 



I. INTRODUCTION 



The characteristic approach has been used successfully 
to carry out numerical simulations of space-times with 
and without sources [l|, 0, y, 0, [a, [a, 0, 0] ■ A significant 
computational effort is nevertheless necessary to extend 
its range of applicability to the simulation of astrophys- 
ically relevant sources of gravitational radiation, such as 
the black hole - neutron star binary problem, where the 
approach can be most useful. As work on the charac- 
teristic formulation to date illustrates p, [^, [lO, [lH , it is 
clear that most three-dimensional characteristic simula- 
tions, even vacuum simulations |9[ are resolution-limited. 
This is particularly true of three-dimensional simulations 
of systems containing compact matter sources [lO| , even 
when the matter source is an extended one. In general, 
all these simulations have been limited in resolution pri- 
marily because of the time required to integrate the equa- 
tions numerically. For instance, at the finest resolution 



'Electronic address: 
1 Electronic address: 
•f Electronic address: 



gomez@psc.edu| 



wbarre to@ula.ve1 



simo@mayu.physics.duq.edu| 



simulation considered in [10i|, tracking a neutron star in 
(close) orbit around a black hole requires approximately 
1.5 months even on one of the fastest processors cur- 
rently available. This is so even though the grid in ques- 
tion (81 X 81 angular points, 123 points radially) is fairly 
moderate by today's standards. To a lesser extent, char- 
acteristic simulations are also limited because of mem- 
ory requirements, although the characteristic scheme is 
particularly economic in this regard. Even though it is 
feasible to equip a single-processor workstation with the 
1.4 Gbytes of memory required by that moderate grid 
size, the time required for the numerical solution on even 
the fastest processor would make such serial simulations 
highly impractical. 

Most of the past code development in the character- 
istic approach [3|, |6[ has been geared towards vector or 
single processor machines. The computational platforms 
available today require instead a parallel programming 
approach in order to perform large resolution simulations 
in a reasonable time, thus a parallel version of the char- 
acteristic code is clearly needed. In the present work we 
show how, with a well thought out yet modest program- 
ming effort, it is not only possible to produce an efficient, 
highly scalable parallel implementation of characteristic 
codes, but to do so in such a way that it becomes straight- 



forward to extend our parallel implementation to new 
physical models. 

We aim for our numerical implementation to be partic- 
ularly useful for long-time simulations of sources of as- 
trophysical interest, which are very demanding in terms 
of the number of grid points on which to advance the 
solution, and thus on the number of floating point oper- 
ations required. With that end in mind, our numerical 
code must scale well in platforms with a large number of 
processors. We show here that our implementation meets 
this goal, making it a potentially valuable tool when ap- 
plied to some of the most interesting astrophysical ap- 
plications of numerical relativity, such as the study of 
black hole-neutron star binary systems, in the close or- 
bit regime up to the tidal disruption of the companion 
star. 

The first astrophysical application we have in mind is 
a characteristic simulation of boson stars in orbit about 
a black hole. For that purpose we calibrate the current 
code for a massless scalar field minimally coupled with 
gravitation. This let us, beyond numerical tests, compare 
and calibrate our code with linear versions of it, radial 
codes, and analytic (perturbative) results reported in the 
literature about quasi-normal modes. The extension to a 
black hole-boson star system is straightforward but not 
trivial, deserving a detailed study of its own and enor- 
mous computational resources. We want to stress that 
the production of gravitational waves by the scattering 
of scalar waves has many mathematical features in com- 
mon with the production of gravitational waves by the 
motion of fluid bodies. 

The article is organized as follows: In Sec. [IT] we re- 
view briefly the standard numerical implementation of 
the eth approach [l|] , and discuss some of its drawbacks 
when applied to high resolution simulations in the char- 
acteristic approach to numerical relativity. In Sec. IIIII 
we present an implementation of the 9 and 9 operators 
which combines their standard description in terms of 
stereographic coordinates with their numerical represen- 
tation on non-conformal grid coordinates on the sphere. 
As this approach differs significantly from our previous 
work, we provide further motivation for this departure. 
Sec. IIVI provides a detailed description of the numerical 
implementation of the approach outlined in Sec. IIIII In 
Sec. |V] we illustrate how the parallel, scalable character- 
istic code framework (LEO) that we have developed, can 
be used to implement the model problem of a scalar field 
minimally coupled to gravity in three dimensions. Sec. I VII 
expands on additional numerical considerations specific 
to the hypersurface and evolution equations, and to the 
boundary conditions, and presents convergence tests of 
our numerical implementation. In Sec. IVIH the viabil- 
ity of the approach is demonstrated by clearly resolving 
several problems which could not be tackled previously. 
We close in Sec. IVIIII with concluding remarks and an 
outline of future work. 



II. THE STANDARD ETH APPROACH IN 
CHARACTERISTIC NUMERICAL RELATIVITY 

The characteristic approach to numerical relativity is 
based on null coordinates a;" = {u,r,x ), with u the 
retarded time, r a luminosity distance and x coordi- 
nates on the sphere [J, i, i, i, i, i, 0, i, [H, [11. In its 
3-dimensional implementation, the angular coordinates 
chosen are stereographic coordinates x^ = {(, () on the 
sphere, which is covered with two stereographic coordi- 
nate patches, as first presented in Ref. [l|] . We summarize 
here the salient aspects of the approach to provide the 
motivation for (and highlight the differences with), the 
implementation described in this article. The standard 
eth approach [l[ is a straightforward numerical implemen- 
tation in stereographic coordinates x^ = {(, C) of the 9, 
9 operators introduced by Newman and Penrose [ij, [l^ . 
Two stereographic patches are used to cover the unit 
sphere, with coordinates (n = tan(0/2)e"^ in the north 
patch and (^s ~ 1/C-/V in the south patch, respectively, 
where (6, (j)) are standard angular coordinates. 

In terms of the dyad q"^ = P{1, i), where P — 1 + CC 
vectors U^ on the sphere are represented by a spin- weight 
1 field U ~ q^Ua (or alternatively, by a spin- weight -1 
field U = q^UA)- This treatment generalizes to tensors 
on the sphere Ta...n+m, which are represented in terms 
of spin- weighted functions obtained by contracting them 
with the dyad q and its complex conjugate q , i.e. 



^ = q^ . . . q " g "+^ . . . q^"+''' 



(1) 



The spin of the resulting scalar function is s = A^ — M. 
Angular derivatives of tensor fields are represented by the 
action of the spin-raising and lowering operators 9 and 
9. For example, the angular derivatives of a vector field 
^ aUb, where ^a are the derivatives compatible with the 
flat sphere metric in the coordinates x , are represented 
by the spin-2 field 9C/ and spin-0 field dU given by 



3C/ - q^q^VAUB, 



W = q^^q^^AUA 



(2) 



The 9 and 9 operators acting on a spin- weight s function 
^ are equivalently defined by 



91' 



P^-'driP-'"^) 



(1 + COac-l- + <*, (3a) 



9* = P^+^dQ (p-'f) = (1 + CO^c^' - sC*, (3b) 
where, in terms of the (real) coordinates {q,p), C ~ 9+*Pi 



dc = do 



idp, di ~ dq + idp. Functions on the sphere 
with spin-weight s transform between patches according 
to 



*/v= (-^ I * 



(4) 



We implement this numerically by laying down a 
two-dimensional grid on each patch, with coordinates 

iqm,Pn), Qm.n = qm+iPn, SUCh that g^ = -1 + (to-3)A, 



Pn = -1 + {n - 3)A, A = 2/{N - 5), and with the grid 
point indices in the range m,n = I . . . N. This grid cov- 
ers the coordinate range —1 — 2 A < {q,p) < 1 + 2 A. 
Ghost zones are used on each side of the grid for the 
discretization of angular derivatives by centered, second- 
order-accurate finite difference stencils. Function values 
at these ghost zones are obtained by interpolation from 
the function values on the opposite patch. In order to 
compute second angular derivatives to second order ac- 
curacy, the interpolations must be evaluated to fourth- 
order accuracy, which can be readily attained using a 
sixteen-point stencil in two dimensions [l|, provided the 
grid covers no less than the range indicated above. For 
some applications Q, we find it necessary to extend the 
grid to a finite overlap, i.e. I?] < 1 -I- e, with e > 2A. 

The set of ghost zones required for the north patch, 
maps onto the south patch (and vice-versa) as per the 
transformation (^^ = 1/Csi va^o a clovcrleaf shape. Re- 
gardless of the number of grid points (or of ghost zones), 
there is a finite overlap between the north and south 
patches. While the points in the overlap area of each 
grid are redundant, we carry them all because we find it 
more efficient to work with rectangular grids. A potential 
problem of doing so is the development of two different 
numerical solutions in the overlap area of each patch, 
only loosely coupled at the stercographic patch edge. 

We now discuss briefly some of the drawbacks of the 
standard "eth" approach when applied to high resolution 
simulations in the characteristic approach to numerical 
relativity, and the motivation for the changes that we 
propose in the next section. 



A. Parallelization of existing characteristic codes 

The first objection that we encounter is in the pro- 
cess of parallelizing our characteristic codes. Because of 
the radial march implicit in the radial integration of the 
hypersurface equations, the natural way to parallelize a 
characteristic simulation is to distribute the angular grid 
among processors, which would integrate the equations 
along a "pencil" of null rays. In the computational "eth" 
approach this means assigning the computation of the so- 
lution over a subset of each stercographic patch to a given 
processor. A similar arrangement, in the context of ax- 
isymmetric simulations, was explored earlier by Bishop 
et. al. [ll. 

Thus, given M x AI processors, we can simply partition 
the N X N stercographic grid on each patch, assigning 
equal square subgrids of extent N/M on each direction 
to each processor. Load-balancing (the requirement that 
all processors in a parallel computation do approximately 
the same amount of work) would in principle be achieved, 
so long as we restricted ourselves to explicit methods, 
thus guaranteeing that the time spent per subgrid re- 
mains constant. The communication pattern imposed 
by the two stercographic patches implementation of the 
"eth" approach [l|] does present an obstacle to effective 



scaling. The mapping of ghost zones at the edge of the 
grid to grid points in the opposing patch is not restricted 
to nearest neighbors. To provide values for these ghost 
zones requires data from a set of grid points whose values 
are scattered among processors in an irregular pattern (in 
the sense that, depending on its location on the angular 
grid, the ghost zones may be obtained from one, two, or 
more processors). This procedure is not only cumber- 
some to program, but intrinsically inefficient. 

If the data required for these ghost zones could be ob- 
tained just from grid points at the edge of an adjoin- 
ing grid, the procedure would simplify considerably. The 
time spent in communication would be substantially re- 
duced and remain constant over the set of processors, 
with a significant impact on the scalability and overall 
efficiency of a code. Unfortunately such an arrangement 
is not possible with a stercographic grid. In addition, as 
we point out at the end of Sec. [Til because of the fixed 
overlap between patches, a significant portion of the grid 
is wasted. While this might be acceptable in small scale 
simulations [0, m, [l3, [U, it needs to be addressed in 
the context of large scale computations as in that case 
it translates into a serious waste of computational re- 
sources. 

Yet another problem that arises from the angular grid 
layout is that of highly non-uniform angular resolution, 
as a direct consequence of using a stercographic grid. 
Considering the expression for the area element in sterc- 
ographic coordinates. 



ds^ 



(1 + CC)^ 



■ dCdC, 



(5) 



it can be seen that there is a marked disparity between 
the resolution at grid points at various places on the 
sphere. Considerable better resolution is attained near 
the equator than at the poles, with the largest disparity 
between a grid zone at the pole (CC — 0) and a point at 
a corner of the grid {(( = 2) where the respective area 
elements have a ratio of 9:1. For some situations, such 
as the case of a matter source in equatorial orbit around 
a black hole, this feature works to our advantage. Con- 
versely, for a matter source in a polar orbit, the matter 
source would be resolved three times as poorly when it 
lies along the z axis (g = p = 0) compared to the res- 
olution obtained when it crosses the equator, which in 
a second-order accurate code translates into a nine-fold 
increase in the intrinsic error in the numerical solution. 
One possible correction for this effect would be to main- 
tain the standard stercographic grid as the basic compu- 
tational grid, but introduce a physical grid related to the 
computational grid by a fish-eye stretch in the angular 
coordinates, and allow the refined portion of the grid to 
follow the compact object. 

A related issue is that of proper resolution of angular 
features, which is critical for characteristic simulations of 
compact objects in orbit around a black hole. In spherical 
coordinates, angular resolution decreases with distance to 
the center, as pointed out in [lO|. There is then a limit 



to the distance at which we can initiahze a characteris- 
tic simulation, a Hniit which is not dictated by physical 
considerations, such as the need to avoid the formation of 
caustics, which would lead to a break down of the coordi- 
nate system. We are also limited by the number of points 
which can in practice be devoted to resolve a companion 
object. 

A proposed extension to the characteristic ap- 
proach [l3| would include an adaptive mesh refinement 
(AMR) strategy, which indirectly would address both 
grid resolution issues mentioned. While this extension 
may prove necessary in simulations of the last stages of 
capture or disruption of a companion star by a black 
hole, the necessary technology has not yet been devel- 
oped in the context of the characteristic approach; to 
our knowledge, applications of AMR in the character- 
istic framework have been made only in simplified one- 
dimensional models [I4I • The intermediate stages of the 
black hole-neutron star problem, where the companion 
remains approximately contained in a finite region could, 
at least in principle, be equally well resolved with fixed 
mesh refinement, an approach successfully used, for ex- 
ample, in [ll,[ll,[23|. The introduction of AMR carries 
with it a whole new set of issues, not the least of which 
is the problem of load-balancing in a massively parallel 
computer. Current characteristic codes present serious 
obstacles to the implementation of an AMR strategy by 
the nature of the angular grid alone. To be effective, an 
AMR implementation would have to be capable of deal- 
ing with refined meshes in multiple coordinate patches. 
To our knowledge, few existing implementations have this 
capability [2l| , and none have been applied in numerical 
relativity. In the present work we present a parallel im- 
plementation on a single distributed grid, and we defer 
the discussion of possible techniques for fixed and adap- 
tive mesh refinement for future work. 



III. AN g OPERATOR BASED ON 
NON-CONFORMAL PROJECTIONS 

A key consideration for the present work is that dif- 
ferent numerical representations of the eth approach can 
be developed by laying down different types of grids on 
the sphere. We consider here an alternative to the stan- 
dard "eth" approach, based on the "cubed sphere" or 
"gnomic" covering of the sphere introduced by Ronchi 
et al. [22| , based on earlier work of Sadourny [2^ . This 
approach is now in common use in global weather simu- 
lations, such as in the MIT General Circulation Model 
(MITgcm), for example [2J]. It has also been ap- 
plied in astrophysical simulations [2^, [2^, in the study 
of wave propagation methods on the sphere [23|, and 
more recentl y, i n the evolution of scalar fields in a fixed 
background [28[. While preparing this manuscript, we 
learned [29J of work being carried out by Bishop et 
al. [30J, on a similar grid arrangement, based on work 
by Thornburg [31 1. 




FIG. 1; The cubed-sphere: covering of the sphere with six 
non-overlapping gnomic patches. 



In the 'cubed sphere' method, a covering of the unit 
sphere with six non-overlapping patches results from pro- 
jecting the sphere from its center onto the six faces of a 
circumscribing cube, whose edges have length two. For 
example, for points on the sphere with Cartesian coordi- 
nates {x,y,z) and angular coordinate — 7r/4 < < n/A 
(where z = cosO), we project the point by tracing a line 
from the center of the sphere though the point {x,y,z) 
to the z = 1 face of the circumscribing cube, determin- 
ing a point with Cartesian coordinates (C/, F, 1), we then 
label the point on the sphere according to the Cartesian 
coordinates {U, V) of its projection on the plane z = 1. 

In order to obtain a covering of the sphere with nearly 
uniform area, we label the points on the sphere by an- 
gular coordinates {a, (i), where U = tan(a), V = tan(/3), 
introducing an equally spaced grid in the angular coor- 
dinates (a,/3), i.e. {ai,(3j) = {iA(;,jAQ), i,j = -N...N, 
A^ = 7r/(47V). Similar projections from the center of the 
sphere to the other faces of the cube provide a covering 
of the sphere with six patches. 

For a finite-difference code, this arrangement is ideal, 
as the grid on each patch is equally spaced in the angu- 
lar coordinates, and the same angular coordinate is used 
on any two patches in the direction perpendicular to a 
boundary. Thus any additional layers of ghost zones in 
the adjacent grid will fall on coordinate lines parallel to 
the boundary. Evaluating function values at those ghost 
zones requires only one-dimensional interpolation along 
coordinate lines parallel to the boundary. For example, 
an A^-th order centered stencil requires A^/2 additional 
layers of ghost cells to be supplied along the edge of each 
spherical cap. The order of the interpolations used to 
supply these points can be selected so as to preserve the 



accuracy and maintain the desired dissipation proper- 
ties of the numerical scheme. In the standard configura- 
tion [22], the spherical caps share common points along 
the edges of the grid, where two grids abut, and at the 
corners of each patch, where three grids meet. These 
common points must have a unique value on each of the 
grids that share them. An approach advocated in the 
literature [23] is to replace the function values at these 
shared points by some form of weighted average. Here we 
dispense with this procedure, by carefully selecting the 
range of the angular coordinates in each patch, in a way 
that precludes the existence of points common to two (or 
three) cubed-sphere patches. 

An advantage of a gnomic decomposition of the sphere, 
crucial to an efficient parallel implementation of the nu- 
merical 3 approach; is that the patches can be laid out so 
that they are non-overlapping. In contrast, in a stereo- 
graphic covering, there is a finite overlap, which does not 
shrink in size as we increase the grid resolution. Although 
it is possible to construct a six-patch stereographic cov- 
ering of the sphere, in which the amount of overlap is 
reduced with respect to the two-patch covering, the area 
of the overlap zone remains constant, regardless of grid 
size. 



Ui = tan(Q:i), Vi = tan(/?i) with -7r/4 < ai./Si < 7r/4). 
The coordinate lines at Ui = const. (Vi = const.) are 
great circles on the sphere which pass through the points 
where the Cartesian axis intersect the sphere. For in- 
stance, the lines at t/3 = const, are great circles spun 
around the x axis, while those of V3 ~ const, are great 
circles rotated around the y axis, with as, (3^ the respec- 
tive rotation angles. Similarly, on each patch, we define 
stereographic coordinates Qi = Ui + ivi, where the an- 
gular coordinates {ui,Vi) on each patch are related to 
Cartesian coordinates by 



{xi,yi,zi) = — (2-Pi,2ui,2t;i), 

{x2,y2,Z2) = — (-2li2,2-P2,2u2), 
^2 

{X3,y3,z3) = — (-2 + P3,-2u3,2v3), 
-fa 

{x4,y4,Z4) = -—{2u4,-2 + P4,2v4), 
^4 

{x5,y5,Z5) = ^(-2w5,2tt5,2-P5), 
^5 



{xe,ye,ze) 



Pfi 



(2ve,2ue 



Pe), 



(7a) 
(7b) 
(7c) 
(7d) 
(7e) 
(7f) 



A. Non-conformal projections of the sphere 

Non-conformal projections are based on non- 
orthogonal coordinates, thus there is no symmetry 
in some operators such as the Laplacian on the sphere. 
It is still straightforward to couple a gnomic grid layout 
with the existing S approach, i.e. while continuing to 
express the 9 and § operators of Ref. [l| on stereographic 
coordinates. In the following we detail how this is im- 
plemented by expressing the angular (stereographic) 
derivatives in terms of angular derivatives in gnomic 
coordinates. A gnomic covering of the sphere is given by 
six coordinate patches 

ixi,yi,zi) = -—{l,Ui,Vi), (6a) 

{X2,V2.Z2) = -^(C/2,1,V^2), (6b) 

{x3,y3,Z3) = 7^ (-1,-^/3,^3), (6c) 

(x4, 2/4,24) = -— {1/4,-1, V4), (6d) 

L)4 

ix5,y5,Z5) = -— {-¥5,1/5,1), (6e) 

{xe,ye,ze) = -— {Ve,Ue,-l) , (6f) 

where {xi,yi, Zi), i = 1...6 are the Cartesian coordi- 
nates of the points, (t/,;, Vi) are coordinates on the sphere 
in the range - \/2 - \ < Ui,Vi < \/2 - 1, and A = 
y'l -I- Vi -|- V^ . (Numerical grid points can be set equally 
spaced in the coordinates (a,/3), related to {JJi,Vi) by 



with Pi = 1 



The gnomic coordinates {U, V) 



and stereographic coordinates {u, v) on each patch are 
related by 



U 



2u 



1 — v?' — v'^' 



V 



2v 



\ — V? — ir 



(8) 



or, in more compact form, in terms of a complex gnomic 
coordinate £^ = U + iV , 



e 



2C 



1-CC 



c 



i 



i + vi + e? 



(9) 



B. The eth operator on the cubed sphere 

We can express the angular derivatives in stereographic 
coordinates, 9^, d? that enter into the 3 and 3 operators 
in terms of angular derivatives {u, v) through the rela- 
tions 



_9__ 

9C 2 \du dv 

and with the Jacobian 



du 



d 



d_ 



d_ .d_ 
du dv 



(10) 



/ dU d£\ 

du dv 

dV dV 

du dv 



1)2)2 



2uv 



2uv 



(11) 



we cast derivatives with respect to stereographic (w, v) in 
terms of derivatives with respect to gnomic coordinates 
(U,V), which are in turn related to derivatives with re- 
spect to gnomic coordinates (a,/3), by 



_d_ 

au 



d 



{l + U^)da 



_d_ 
dV 



d 



{1 + V^)d(3' 



(12) 



The spin-raising and lowering operators 9 and 3 |32l | act- 
ing on a spin s function ^P, Eq. ([3]), can be written as 



g* = (i + cc) 

+ <*, 
S* = (i + CC) 



1 a* i d-^ 



1 + C^ 3a 1-C d(3 



1 d^ 



a* 



l + C da 1 - (2 dp 



(13) 



(14) 



where the values of (C, C) for each grid point are computed 
from the respective {£,,£,) values as per Eq. ©. 

To complete the prescription of the 3 operator on the 
cubed sphere we must specify the transformation rule for 
spin-weighted functions on the sphere. The stereographic 
coordinates on the various patches transform according 
to 
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(15a) 



(15b) 



(15c) 



(15d) 



where Eqs. (J15a[) relate neighboring equatorial patches, 
Eqs. (J15bp and (|15cp supply the coordinate transforma- 
tions between equatorial patches and the north and south 
patch, respectively. Eqs. (|15dp . which relate diametri- 
cally opposed patches on the sphere, are not strictly nec- 
essary for a numerical implementation and can be de- 
duced from (|15bp and p5cp . Our convention for the 
gnomic parameterization follows Ref. [22| . and it has 
the advantage that the orientation of the (u, v) axes on 
each patch is chosen so as to reduce the amount of book 
keeping needed to transfer information between patches. 
From Eq. (|15p . adjacent patches with coordinates Ci and 



Q are related by AQ = {BQ - l)/iBQ + 1), where 
A = l,±i and B = 1,+*, while opposite patches are 
related by Q — C/Cj, with C = ±1, in particular, we 
recover the coordinates transformation between "north" 
(Cn = Csj and "south" (Cs = Ce) patches, Cn = 1/Cs, as 
in Ref. ^. 

Given the dyad q^ = P{l^i), P{CX) = 1 + CC, and 
expressing q = q^dA = P{C, Qd^ in two coordinate sys- 
tems x^ = (C, C) and x^ = (C'lC'); i* can be seen that 



where 



wiCO 



q = W{C,C)q', 

p(C,C) dC 



(16) 



P(C',C')5C' 

and where the substitution (' — C'(C) is understood in 
the right hand side of (fT6|) . From this we deduce the 
transformation rule for spin-s functions on the sphere, 
given here only for the case of adjacent patches 
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For the case of functions with spin-weight zero the trans- 
formations reduce to "^jXjXj) = ^i(CijCi)- 



IV. NUMERICAL IMPLEMENTATION 

We give here a summary of the numerical techniques 
used so far in the LEO framework. It is worth noting that 
the framework is easily extensible, and thus we are not 
restricted, for instance, to the particular choice of radial 
grid made here, nor to the choice of radial or time integra- 
tion schemes used in the present work. In the subsection 
on finite-difference operators on the sphere, for instance, 
we describe higher-order extensions that we have elected 
not to use in the example application considered here, as 
they are inconsistent with the radial and time integration 
schemes, which we have taken unchanged from [3|. 



A. Radial grid and finite difference operators 

Following [3^, we take the computational radial grid 
to be equally-spaced in the compactified coordinate x = 



r/{R + r), restricted to the range a;_B < a^ < 1, i-C. Xk = 
XB + {k- l)Ax, k = l,...,N^,Ax={l- xb)/{N^ - 1), 
with xb ~ TB/iR + rB)- We express radial derivatives 
in terms of the conipactified grid Xi, via the relation 
dx/dr = (1 — x)'^/R, e.g. 






(l-^fe+^)'(/fe+l-./fc) 



R 



Aa 



(l-Xk)^ {fk+i - fk-l) 



R 



2Ax 



(18) 
(19) 



B. 



Centered finite difference operators on the 
sphere 



We construct an equally-spaced grid on the gnomic 



coordinates x = {a, (3), with 
/3, = -V4+(j-i)A, andA = 



N^. 



= -7r/4 +{i- i)A 

7T/{2N^),t,j^l 

The useful part, exclusive of ghost zones, for each of the 
coordinates ranges from — 7r/4 + A/2 to 7r/4 — A/2. With 
this arrangement, the points with |u| = 7r/4 or \v\ ~ 7r/4 
are excluded, and thus we avoid storing double values for 
the points at the edges of each patch, and triple values 
for the points on the corners where three patches meet. 
Adding A'^ ghost zones on each side of the grid allows us 
to evaluate derivatives to order N 
stencils of the form 



2Ng with centered 
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FIG. 2: Convergence rate of the 9 operator, built upon angu- 
lar derivatives of order 2, 4, 6, and 8 (indicated in the graph 
as circles, squares, diamonds and triangles, respectively), act- 
ing on 2Y43, and with grid sizes ranging from A*',; = 16 to 
Nc = 128. 



C. One-dimensional interpolation of ghost zones 
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^^Ck{f, 
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fc=i 



The coefficients for the derivatives, up to 8-th order, are 
given in Table [D We can verify that the coefRcients of 
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1/2 
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8/12 


-1/12 
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3/4 


-3/20 


1/60 
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4/5 


-1/5 


4/105 


-1/280 



TABLE I: Coefficients for centered angular derivatives. 

Table |I] for each of the derivatives are correct by not- 
ing that the numerical error of derivatives of order N is 
within the level of round-off when applied to a polyno- 
mial test function F of order N or lower. We have also 
verified the proper convergence rate of the numerical 3 
and 9 operators when applied to spin-weighted spheri- 
cal harmonics [9[ . Fig. [5] shows the proper convergence 
rates of the 3 operators constructed from derivatives of 
second, fourth, sixth, and eighth order when applied to 
the spin-2 spherical harmonic 2^43 on grid sizes ranging 
from Nc = 16 to Nc = 128. 



Since the ghost zones required to evaluate derivatives 
fall on coordinate lines parallel to the boundary, we can 
obtain function values at these ghost zones with one- 
dimensional interpolations. We use standard Lagrangian 
interpolation formulae to A'^-th order accuracy. 
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j¥-' 



(Xi Xj) 



(21) 



adapted to equally spaced grids, i.e. Xj ~ xq+j A. Fig. [3] 
shows the calibration of the interpolation routines with a 
test function consisting of a polynomial of order 15, i.e. 



N 

PN{a) = y^ Cj g' 



i=0 



< a < - 

4 - - 4 



(22) 



H 



C 



for A^ = 15, where the coefficients Ci,i = 1 . . . N are cho- 
sen randomly, subject to the condition \ci\ < 1. The in- 
terpolants display convergence to the correct order (3,5, 
7, and 9-th order, respectively), for grid sizes in the range 
8 < A^C ^ 256. For this range of values (8 < A^; < 256), 
there are from 32 to 1024 points in the great circles deter- 
mined by the intersection of the sphere with any of the 
Cartesian coordinate planes. As expected, for smooth 
data such as our test function, for sufficiently large grid 
sizes, the error goes down to round-off level when using 
the higher-order schemes. This saturation effect is al- 
ready visible in the plot for the 9-th order interpolator 
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order accuracy by evaluating the area element in gnomic 
coordinates, 






(23) 



evaluating the function value on grid cell centers, ftj 
f{Ui,Vj), and summing over grid cells, 



TV,. N, 



r ^ ^ 1 



^^, (l + t/2)(l + V;2) 



tU''''{i + uf + vj^)y^ 



A^ 



(24) 



where A stands for the grid spacing on both coordinates 
(a,/3), which we have taken to be the same. Note that 
since the spherical patches do not overlap, the integral 
over the sphere is just the sum of the integrals over the 
individual patches. Fig. |4] shows the converge of the in- 
tegral of the area element itself to the correct answer 
oi J dfl = 47r for grid sizes in the range of A^^ = 8 to 
A^^ = 512. The measured convergence rate is 2.0, in full 



FIG. 3: Convergence rate of the various interpolation schemes 
used. Shown in the graph are the 3'^'' order (circles), 5"' or- 
der (squares), 7"* order (diamonds), and 9*'' order (triangles) 
interpolators. For the highest order interpolator used (9-th 
order), the error goes down to double-precision round-off level 
(~ 10^^^) when more than A'^^ = 100 angular points per patch 
are used. 



when angular grid sizes reach approximately N,^ = 100. 
As indicated in Sec. Illli we have chosen the range of the 
spherical coordinates {a, f3) so that there are no over- 
lapping points at the edge of each patch. This avoids 
the awkward procedure of averaging values from differ- 
ent patches to obtain a single- valued function throughout 
the computational grid. 

The number of ghost zones, the order of the finite dif- 
ference approximations and the order of interpolation, 
while related, are not directly tied to each other. One re- 
quirement is that we must have enough ghost zones (Ng) 
to compute the finite-difference approximation to the de- 
sired order, Np; and since in general we want to use 
centered differences, the relation Np < 2Ng must hold. 
If we wish to maintain the symmetry of the interpolation 
stencils, N] < 2Ng + 1 must also hold. For the cases we 
have considered, we find that our algorithms are stable 
if Nj > Np + 1, with the inequality required only in the 
case of Np = 2, the lowest order of finite-differences that 
we considered. We are otherwise free to vary the number 
of ghost zones as dictated by efficiency considerations. 



D. Integrals over the sphere and volume integrals 

Integrals over the sphere and volume integrals arise 
naturally, in particular when computing norms of various 
quantities. We evaluate integrals on the sphere to second 
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FIG. 4: Convergence rate of the integral of the area element 
over the sphere, for grid sizes ranging from A^^ = 8 to N^ = 
512. The markers indicate the error of the area element at 
the corresponding resolution, the line is the least-squares fit, 
yielding a convergence rate of 2.0. 

agreement with the expected result. 

Volume integrals are computed similarly to integrals 
over the sphere, but in this case evaluating the spherical 
contributions mid-point in between radial points, i.e. 



dn 



f^'dr=^^- 



^^(l + [/f)(l + ^2)^^ 



tf^,i^ + uf + vf)y^ 
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Ax. (25) 



We replace the flat volume element, dV = r'^drdH., with 
the volume elem.ent corresponding to a Bondi metric, 
dV = r^e'^^dr d^, when appropriate. To speed up the 
evaluation of integrals, we pre-compute the area element 
on the sphere, Eq. ((23)) . 



E. 



s function, 



Accuracy of the spin-weighted spherical 
harmonic decomposition 



We make use of spin-weighted spherical harmonics 
sYirn throughout this paper, following the convention 
of [9| . In order to estimate the error introduced when we 
perform a spin-weighted spherical harmonic decomposi- 
tion, we look at how well the orthonormality condition 



Ylm. ..Ylmd^ = 5l_v5, 



(26) 



is preserved (at the numerical level) for spherical har- 
monics with spin- weight s = 0, 1 and 2, for a range of 
values of £ and m and angular grid sizes. As expected, 
the numerical value of the integral converges to the ana- 
lytic result to second order on the grid spacing, since we 
have chosen to use a second-order integration algorithm. 
Fig.[5]illustratcs one instance, where we have taken s ~ 0, 
£ =: 6, with m = —6 ... 6, and varied the angular grid size 

We can also place an esti- 
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FIG. 5; Convergence of the orthonormality condition, 
illustrated here by computing the convergence rate of 
/goyimoyim = 1, for the case / = 6, m = . . . 6, on grid 
sizes ranging from A'^^ = 32 to A'^^ = 64 



mate on the accuracy of the projection of a spin-weight 



Clm{F\ 



F sYlrndn 



(27) 



based on the magnitude of the off-diagonal values in 
for a given grid size. When projecting the test functions 
Yi' ra' into the spherical harmonics Y; „ for Z = . . . Imax-, 
m = —I . . .1, ai the analytic level we would expect to ob- 
tain zero for all coefhcients, except for c/' m' which would 
be identically one. We find that grid sizes of Nq = 64 and 
larger are sufficient to keep the error in the coefficients 
to within one part in 10^, which again is consistent with 
our integration scheme being second-order in the angu- 
lar discretization. Fig. [6] shows the error in the coeffi- 
cients computed for yi'm/, Z' = 6, m' = 3 on a grid with 
N(; = 64 points. We have omitted from the graph those 
coefficients for which the error is already at the level of 
round-off. 
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FIG. 6: Error in the coefficients c; „ when the function being 
projected is the spherical harmonic lea, on a grid of A'^^ = 64 
points. Coefficients whose error is at round-off level are not 
shown. 



The preceding description of the numerical implemen- 
tation is complete but for one key aspect, namely our 
parallclization strategy. In our framework, the six cubcd- 
sphcrc grid patches arc decomposed into computational 
sub-patches, each with the same number of points on 
the angular directions, for efficiency reasons. These 
sub-patches are distributed among processors, and the 
ghost-zone values required for the computation of angu- 
lar derivatives are communicated by the use of message- 
passing calls [3J| . The radial direction is not distributed, 
as the characteristic algorithm requires a radial march 
for the integration of the hypersurface equations as well 
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as the evolution equations [33[ . The computational com- 
plexity of a parallel implementation via message pass- 
ing lies in that, knowing the location of its assigned grid 
sub-patch on the global grid, each processor must deter- 
mine which processors are its nearest neighbors, i.e. to 
which processes it must supply ghost-zone information 
(and also receive that information from) . Due to the rel- 
ative orientation of the cubcd-sphcrc patches, we need 
to know whether the order in which the ghost zones are 
traversed must be reversed for sub-patches on the edge 
of a cubed-sphere patch. Since the sub-patch to proces- 
sor mapping remains constant during a simulation, the 
relevant information needs to be computed only once, 
and at any rate, it incurs no measurable overhead in the 
computation involved in a simulation. An efRcient and 
scaling implementation of the message passing itself re- 
quires only a small subset of the full MPI functionality: 
a few calls to set up the appropriate groups of processors; 
sends, receives and waits (for ghost zone communication); 
some additional reduction operations (to accumulate in- 
tegrated values), and some broadcasts (to propagate pa- 
rameters). Exclusive of file access operations, only 14 
MPI functions in all are invoked. 

Having established that all the key computational as- 
pects of the framework are in place, and have been cor- 
rectly implemented, wc proceed next to consider specific 
applications of the framework to systems of physical in- 
terest. 



V. A THREE-DIMENSIONAL MASSLESS 

SCALAR FIELD SCATTERED OFF A 

SCHWARZSCHILD BLACK HOLE 

We use the numerical formalism developed in the pre- 
ceding sections to solve numerically a model problem con- 
sisting of a self-gravitating massless scalar field in three 
dimensions. Our starting points are Ref. [SJ] for a descrip- 
tion of the vacuum problem, and Ref. |ll| for the coupling 
of the scalar field to the gravitational metric fields. We 
use coordinates based upon a family of outgoing null hy- 
persurfaces, and we let u label these hypersurfaccs, x^ 
{A = 2, 3) label the null rays and r be a surface area 
coordinate. In the resulting x" = {u, r, x ) coordinates, 
the metric takes the Bondi-Sachs form [Sa, [3611 



function J = hAsq^q^ /2, and by the real function 
K = hABQ^q^ f^, where K"^ = 1 + J J. The metric func- 
tions U^ are similarly encoded in the complex function 
U = U^QA- Thus, it is necessary to introduce the inter- 
mediate spin- weighted variable Q ~ QaQ , as well as the 
(complex differential) operators 3 and 9 (see [1| for full 
details). 

Treating the Einstein-Klein-Gordon model problem 
consistently within the LEO framework requires some 
modifications to |ll| , specifically to the wave equation for 
the scalar field (Dt/) = 0) which is given by Eqs. (21)-(27) 
of [11|. We substitute all second-order angular deriva- 
tives of the metric fields in terms of 9 and 9 operators 
acting on the additional fields i^ — 9J, k = dK and 
B = (5(3 introduced in Ref. Q, whenever possible. A 
consistent treatment is obtained by introducing the ad- 
ditional variable 



V' = 3x, 



(29) 



where x — ^^j so that the scalar field equation is also in 
first-order differential form in the angular variables, on 
a par with the approach of [8| for the metric equations. 
The Bondi-Sachs hierarchy of hypcrsurface equations. 



1^,7 



B.r 



r^U, 



(r^W), 



9Jr, 



(j.rJ.r - Kl) + 2Trr{(t>^rf 



= 3Ar, 



(30) 
(31) 
(32) 

(33) 



— K{k_r + l^,r) + ^J.r + J^J,r + ^K^r 



^^ [P {J,r ~ J^J,r) + 9J (J. - JV,.)] 

+2r^B.r - 4rB + 167rr0,rV , (34) 

e^P {KQ - JQ) , (35) 

5r| e2/3 C^-K (SB + BB)+J (9B + B'-) 



ds^ 



e^P{l + W/r) - r'hABU''U"]du' 



2r^hABU"dudx' 



f liABdx dx , 



2e^f^dudr 
(28) 



where W is related to the more usual Bondi-Sachs vari- 
able V hy V = r + W, and where h^^hBc — ^c ^^'^ 
det{hAB) = det{qAB), with qAB a unit sphere metric, 
given in terms of a complex dyad qa satisfying q^qA = 0, 
q'^qA = 2, q^ = g^^gs, with q^^qBC = 5^ and 



^c 



qAB = \{qAqB + qAqB)- 

ate variable Qa = r^e^'^'~ hab'- 

sors on the sphere by spin-wighted variables [l|. The 

conformal metric Hab, is represented by the complex 



We also use the intermedi- 
U^. We represent ten- 



^,r 



{i^ - k)B) -1 + 2 rW + ySt^.r 
r* _ _ 1 

-e-^'^ — U^r{KU,r + JU,r) > 

27r^- [2KiJj'iP - Ji? - J-i}?] , 



(36) 

(37) 



now includes an additional consistency condition, 
Eq. (P?)) . and the equations for /?, Q and W ~ W/r"^ 
are modified to include the source terms as shown above. 
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The evolution equation for the metric field J is given by 

2 irJ),ur ~{r~^V (rJ)^^^ = -K {rdU^r + 2dU) 

+ -e2'3 (35 + B^) - (^rWr + W^J+Jh + JPu 

+ ^e2V, (38) 

with the quantities TZ, Jh and Pu as in Eqs. (24)-(26) 
of [8|. The scalar field evohition equation follows from 
Eq. (21)of[Il|, 



V 

^X,ur I X.7 

r 



The source term N^ is 

= 2/3 r I 



X 



+ N.. 



(39) 



iV^ = — 



r 



— ( jgv + J^i') + — SV' 

2r ^ o. 



r 



KB - JB 



^{KQ-JQ)-'^ 



+ i^(^"- + ^^))? 



KB-JB-^{KQ-JQ)-'^ 



+i^(^^+^»)?; 



-i(t/V + Ut/A - -(pJSU + dU) 
r 2 ' 



(40) 



Following [iJl, we have used the shorthand /i = 9 J, 
and eliminated the radial derivatives C/,r and U^r using 
Eq. dSll), 



Q^r^e~''''{KU,r + JU,r). 



(41) 



The data required on the initial null cone are the evo- 
lution variables J and (j)- Given boundary values at a 
fixed value of r, the remaining variables {v, k, /3, B, Q, 
U and W) can be determined on the initial null cone 
by explicit integration of the hypersurface equations (see 
[m for details). The evolution equations ([55|l and ([55]) 
can then be used to find J and 4> on the next null cone, 
and the process repeated to determine the spacctime to 
the future of the initial slice. 



Eqs. (|37|) and ([39)) . For a Scharwzschild background, 
the metric fields J, /3, U, v, k and B are zero, and 
V ~ r — 2M. The source term in Eq. (|40p reduces to 
N^ = dip/r, and we are left with the system 



2X,. 



2M\ 



2Mx^dl 



i^,r = (^X; 



(42) 



For the simulations we discuss in the present work, we 
will be interested in solutions of the scalar field on a 
fixed background with definite angular dependence, as 
discussed in the next sub-section. 



B. Quasi- normal modes in a Schwarszchild 
background 

The linear equation for the scalar field on a fixed back- 
ground, Eq. (|42|) is separable, i.e. its solutions can be 
written in the form 



hr,x ) = 2^ 2^ Xim[u,r) 



(43) 



£=0 m=-e 



with the X coordinates in the sphere, and where each of 
the Xim satisfies the one-dimensional wave equation in 
the plane (u, r) 



2X, 



1 



2M\ 



X.r 



r ) 

2M (.{i^X) 



(44) 



where we have used the property 9§x = — ^(^ + l)x pTj . 
Eq. (|44|) is the usual equation governing the scalar per- 
turbations of a Schwarzschild black hole [38| , written here 
in characteristic coordinates (w, r, x ). It can be put in a 
more familiar form by writing it in the coordinates (t, r*), 
with u = t — r^:, and where r» is the usual "tortoise" co- 
ordinate, r^ =r + 2M\n{r/2M ~ 1). 

X,u - X.r.r. + V{r)x = 0, (45) 

where x = ''^ and the potential V{r) is given by 

iie + 1) 



V{r) = 1 



^) 



2A/ 



X, 



(46) 



A. Scalar field on a fixed background 

The above system of equations describes a self- 
gravitating scalar field. In the limit of small amplitudes, 
\(j)\ « 1, the scalar field can be treated as a pertur- 
bation propagating on a fixed background. This consid- 
erably simpler model is contained in the fully nonlinear 
case, and is implemented in our code by integrating only 



and we have denoted it by V to avoid confusion with 
Bondi's V which we use throughout this paper. Eq. ([45]) 
has been studied extensively ^38., :39;i40], its most salient 
feature being the existence of quasi-normal modes, whose 
frequencies have been tabulated, see for example [4(j |. 
Note that the right-hand side of Eq. (|1^ is the correct 
form of the potential in (u, r) coordinates. It differs by 
a factor of (1 — 2M/r) from the potential as given in 
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Eq. (|45p . see [38[, because that factor is precisely the 
Jacobian of the coordinate transformation, d^/dr, = 1 — 
2M/r. 

In the remainder of the present work we will use both 
the quasi-normal mode equation, Eq. (j44p , and the linear 
system, Eq. (|42]) , as tests of the validity of our numerical 
implementation. We do this in an incremental fashion, 
solving Eq. (J44|l for fixed values of £, and comparing the 
effectiveness of the numerical integration scheme and of 
our boundary conditions in reproducing the quasi-normal 
modes. To this end, wc implement a purely radial code 
for Eq. ([m that employs the same numerical integration 
scheme that is used in the "linear" code (which solves 
Eq. (gH) and Eq. dH])), and in the full nonlinear code. 
Using this radial code, we can isolate the effects arising 
from the inner-boundary treatment at r = 2M by imple- 
menting Eq. (|44)) as indicated, in outgoing null coordi- 
nates, using both a non-compactified coordinate r, with 
a simple extrapolativc boundary condition at the outer 
boundary Tout > 2 A/, and the compactificd coordinate 
X = r/{r + R), where the outer boundary lies at future 
null infinity. The use of a non-compactified coordinate 
allows us to isolate any effects that may arise due to the 
non-uniform coordinate velocity introduced by the com- 
pactified coordinate x. Conversely, simulations using the 
compactificd coordinate avoid the effects of placing the 
outer boundary at a finite distance. 

We also implement the equivalent of Eq. (|^l)l in ingoing 
null coordinates, {v,r), with f = < + r*, namely 



■^ X,vr 



2M, 

(1 )X. 

r 

2M £{£+!) 



(47) 



In ingoing null coordinates, the slices at u = const pene- 
trate the event horizon r = 2M , effectively providing for 
an excision scheme, where evolution can be stopped at a 
finite number of points inside the boundary, because the 
behavior of the field inside the horizon does not affect the 
solution outside. Evolutions in ingoing coordinates are 
carried out on a non-compactified radial grid, for which 
boundary data are required at a fixed value of rout > 2M. 
Because of the presence of this outer boundary, simula- 
tions in ingoing coordinates can only be run for a limited 
time, typically u ~ 2 rout, before outer boundary effects 
influence the signal extracted. A similar effect is seen 
when using outgoing, non-compactified null coordinates. 
When using compactificd coordinates, no such effects are 
seen, as expected. A detailed comparison between in- 
going and outgoing versions of characteristic systems of 
equations, in compactificd as well as non-compactified 
coordinates, along with their relative advantages and dis- 
advantages for specific applications, is worthwhile but lies 
outside the scope of the present work and will be reported 
elsewhere. We will refer only briefly to these issues in the 
remainder of this work. 



C. Energy carried out by the scalar field 

As a useful physical indicator we calculate the balance 
of the scalar field energy contained between the inner 
boundary and null infinity. The expressions we give here 
are valid in the linear case, where the background metric 
is that of Schwarzschild. For a more general approach 
to this issue, the linkage integrals have to be calculated, 
specifically the asymptotic Killing vector field must be 
parallely propagated from null infinity j4l| . 

Restricted to the background case then, given a Killing 
vector field ^'^ of the metric g^i,, £(_g^ij — 0, we can define 
the conserved quantity 



c^JnTd^, 



(48) 



In particular, selecting the time-like Killing vector ^'^ ^ 
d'^u, and for a surface of constant u, C is the energy con- 
tained on the surface. 



E{u) 



T^dV, 



(49) 



where dV is the volume element of the surface at constant 
u. For a sphere at constant r, C represents the energy 
fiux across the surface, 



p{u)^ frrydn, 



(50) 



with dfl the solid angle element. The relevant compo- 
nents of the stress-energy tensor for a massless scalar 
field are 



rpU 



1 
4^2 



V_ 
2^ 



j(g<? 



o' + jfa 



1 



-2/3 



<j>.r (C/g0 + Ud(t)) 



rpr 



-2/3, 



cb,u~j(l),r + ^{Udq^+UB(l)) 



(51) 



(52) 



In the case of a linear scalar perturbation on a 
Schwarzschild background, the energy content of a hy- 
persurface at constant u is given by 



E{u)^\J (^1 - ^^ (r0,,)' + S^S? 



drdn . (53) 



The power radiated at time u across a surface of constant 
r, such as the inner boundary, which in our simulations 
we place close enough to the Schwarzschild black hole, is 



Pin{u) 



'-^) 



r^dn 



(54) 



For the flux across the inner boundary, the integral as 
well as the spatial and time derivatives are to be taken 



13 



as evaluated at r = ri„ . Analogously, the power radiated 
at time u at null infinity is the limiting form (as r — > oo) 
of the above expression, i.e. 



Poutiu) = / {r<j),u)^dn 



(55) 



where we have used the behavior of the scalar field near 
null infinity J^ to simplify the expression. With these 
definitions, the following global energy conservation law 
holds 

T.{u)^E{u)+ / [Pout{u') - Piniu')]du' = const. (56) 



Even though the expressions given above hold only in 
the limit in which dt is a Killing vector of the metric, 
we expect them to hold in an approximate sense for our 
nonlinear evolutions, so we use them as a criterion for 
code testing. 

As stated previously, we use the radial code to calibrate 
the fully nonlinear, three-dimensional LEO code in the 
linear regime. When computing the energy in the radial 
code, we make use of the property 



iQchdfl 



>Bdd>dD.. 



(57) 



(see [33 ) ■ Since the data we pose are pure spherical har- 
monics, the integral in the right-hand side is proportional 
to the norm J (j>(f> dfl. Eq. ((57|) allows us then to properly 
account for the contribution of the angular derivatives of 
the field to the energy (|49|) when using only the radial 
code. 



VI. ADDITIONAL NUMERICAL 
CONSIDERATIONS 

A. Hypersurface equations 

The integration of the hypersurface equations does not 
present any inherent difficulty as they are discretized at 
mid-point between grid points as per [8|, |ll| . An impor- 
tant issue which arises because of the parallel implemen- 
tation of our algorithm is that after each step in the ra- 
dial march, that is, after each hypersurface equation has 
been advanced radially one grid point, we must synchro- 
nize the variable which has just been integrated. By this 
we mean that we communicate the ghost zone values to 
the processors carrying out the integration in neighboring 
patches. Since communication is an expensive operation 
even on the most tightly coupled parallel computers, we 
take the approach of explicitly synchronizing a variable 
only if an 9 (or 9) operator will be applied to the vari- 
able in question. An alternative approach would be to 
incorporate the synchronization into the 9 (and 9) op- 
erators. The first approach requires more book-keeping 
on our part, whereas the second is more straightforward. 
Because of the number of 9 (or 9) operations that appear 



in the full nonlinear equations, however, the performance 
difference between these two approaches is significant. 
For this reason we take the first approach, reducing to 
the minimum possible the amount of communications, 
with a substantial increase in performance. 



B. Evolution equations 



The evolution equation ((38|) for J is treated as reported 
in [ll| . except that the first two radial points are subject 
to the boundary condition explained below. The evolu- 
tion equation for the scalar field is recast in terms of the 
two-dimensional wave operator 



where x = ^0 and Eq. p9|) reduces then to 



where 



n = -{W/r),rX/r + N^. 



(58) 



(59) 



(60) 



Since all two-dimensional wave operators are conformally 
flat, with conformal-weight —2, we can apply to (|55|) a 
flat-space identity relating the values of x at the four cor- 
ners P, Q, R and 5* of a null parallelogram A, with sides 
formed by incoming and outgoing radial characteristics. 
In terms of %, this relation leads to an integral form of 
the evolution equation for the scalar field 



Xq=Xp + Xs- Xr + 



1 



dudrJi. 



(61) 



The corners of the null parallelogram cannot be chosen 
to lie exactly on radial grid points, thus the values of x 
at the vertices of the parallelogram are approximated to 
second-order accuracy by linear interpolations between 
nearest neighbor-grid points on the same outgoing char- 
acteristic. Approximating the integrand by its value at 
the center C of the parallelogram (evaluated using aver- 
age values between the points P and 5), we have then 



XQ = XP + XS- XR 

Au . 
+ —r\rQ-rp + rs- rn) Tic- 



(62) 



The evolution algorithm for the metric function J follows 
the procedure outlined in [3, [a, llH . As with the hyper- 
surface equations, we synchronize the fields (j) and J, i.e. 
we communicate the ghost zone information from each 
patch to their neighbors, immediately after advancing 
radially these two fields with their respective evolution 
equations. 



C. Boundary treatment for the evolved fields 

For the ingoing formulation, we set the field values 
(u,r — rout) = 0, and we march inwards until a few 
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points beyond the black hole horizon (r ~ 2M). Since 
the past light cones tilts outwards once inside the horizon, 
the values computed just inside the horizon can never af- 
fect those points of the grid that lie outside. This scheme 
provides then an extremelv simple and effective form of 
excision, as discussed in 0, I5l .l42l in the context of char- 
acteristic evolution, and in [43| in the context of 3+1 
simulations in the Bondi-Sachs gauge. 

For the outgoing formulation on a non-compactified 
radial grid, we use simple extrapolative boundary condi- 
tions at the outermost point, i.e. the field x at the last 
point is set equal to the value of x ^^ the point imme- 
diately before. This approximation is justified for suffi- 
ciently large r as the field 4> behaves, to leading order, 
as ~ 0{l/r). Our treatment of the inner boundary is 
motivated by physical considerations that arise naturally 
in the study of quasi-normal modes. It can be seen from 
Eq. (|15)) that when the potential V{r) goes to zero, as 
it does in the limits r — > 2M and r ^ oo, the solutions 
to Eq. (|45|) are traveling waves, x- = FL{t + r^,) and 
X+ = Fait — r»). In the linear approximation then, it 
is consistent to apply an open boundary condition to the 
the scalar field (j) based on the assumption that, at the in- 
ner boundary, the field behaves as a left-travelling wave, 
X = FL{t + r^,). It follows also that the same condition 
must be applied to the spin- weight 2 metric field J. In the 
linear approximation, Eq. ([55]) reduces to Eq. ([15|) . with 
the potential V{r) corresponding to that of a spin- weight 
2 field, see [38| . This open boundary condition is equiv- 
alent to stating that the fields x~f4> and rJ propagate 
towards the horizon along the incoming characteristics 
of the two-dimensional wave operator, Eq. (j58p . In prac- 
tice, we implement this condition for the first two points 
of the radial grid, and use the evolution equations for x 
and J elsewhere. 

In the non-linear case, the horizon can no longer be as- 
sumed to be static, rather it is dynamically distorted and 
grows as the scalar field accretes into the black hole. Our 
boundary condition is applied always to the same set of 
points, which are subsequently enveloped by the growing 
horizon, thus any inaccuracy we might have introduced 
in those first two points can not have any effect on the 
exterior spacctime. 

Our approach suggests the following iterative method 
to treat the inner boundary, in a manner which is con- 
sistent with the open boundary condition: (1) as a first 
approximation, solve the homogeneous equation (j45p for 
the first two radial points, i.e. assume the evolved fields 
propagate along incoming characteristics up to the re- 
tarded time u -f- Am, and (2) with the values predicted 
for the fields at time u -I- Au, correct the right-hand side 
of the full evolution equations. 



D. Tests of second order convergence 

The simulations for the tests reported in the remain- 
der of this section are conducted in compactified outgoing 



(retarded) null coordinates. To verify that the numerical 
algorithm is globally second-order convergent, we com- 
pute the L2 norm of the relative residuals for three grid 
sizes, e.g. 



= [Xc-Xm]'^dxdn, 



Qmf 



Xn 



XfYdxd^, 



(63) 



where the c, m and / subscripts denote the field as com- 
puted on coarse, medium and fine grids, respectively. 
The field is evolved from an initial retarded time u = 
and the integrals ((63)) are calculated at the same final re- 
tarded time u, using the same set of spatial grid points, 
obtained by appropriately subsampling from the fine and 
medium grids to the coarse grid. Here wc take the angu- 
lar (and radial) grids to be in a proportion of 1 : 3 : 5. 
Grids in these ratios have a common set of points that 
align directly, and thus do not require interpolating cell 
values from the finest to the coarser grids. In this case, 
given the values Qcm and Qmf , it can be shown that the 
order of convergence 0(A") of the algorithm can be read 
by solving for n the following equation 



{Qcm/ Qmf) 



1/2 



1 - 1/3" 
1/3" - 1/5" 



For this test we evolve the initial data 



x(0,r,x^)=Ae 



-(r-rof/a' 



Y,n 



(64) 



(65) 



with A = 10^^, whose the radial profile is character- 
ized by rg = 3M, a = 5M, and whose angular depen- 
dence is given by £ = 4 and to = 2, from u = up 
to u = \M . We perform three simulations, on the an- 
gular grid sizes A^^ ~ 10, 30, 50 and the corresponding 
radial grid sizes N^ = 501, 1501, 2501, for which wc take 
152, 456, 760 time-steps, respectively. From ()64p we find 
that the measured order of convergence is n = 2.05, in 
excellent agreement with the expected second-order con- 
vergence. It should be noted that this procedure tests 
the Cauchy convergence of the code, providing a basic 
check of the consistency of the discretization. For low 
amplitudes (in the perturbative regime), and for a given 
value of £, the scalar field profiles computed with the 
fully three-dimensional code match, to within second or- 
der, the profiles obtained with a purely radial code which 
solves Eq. ([44]), as expected. 

We want to stress that the boundary conditions, the 
initial data and the marching algorithm for the scalar 
field used in this numerical test are all the same as those 
which we have used to calibrate the radial code, the solu- 
tions of which we use here in place of an analytic solution. 
In fact, the convergence rate for the radial code is exactly 
2.00 for the radial grid sizes oi N^ = 501, 1501, 2501, mea- 
sured at u = IM with its respective subsampling, as per 

Eqs.dSSI-dMD- 

For sufficiently low values of (£, tti), the angular grid 
sizes N^ used in the convergence test are adequate. For 
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a given angular grid size, it is also possible to reduce 
the angular error by increasing the order of the angu- 
lar derivatives, for example, to fourth order or higher. 
The increased computational expense is offset by the in- 
creased accuracy obtained; in a parallel application there 
is also the potential for additional overhead because more 
ghost cells must be communicated. In practice we ob- 
serve that, for the smallest angular grid size considered 
{N^ ~ 10), changing the discretization of the angular 
derivatives from second to fourth-order increases the ex- 
ecution time by about 20%. 

In the work reported here, since the radial and time 
integration are carried out with a scheme that is second- 
order convergent, we have opted to use second order accu- 
rate angular derivatives, as with this choice the nonlinear 
code exhibits second-order accurate Cauchy convergence. 

For the initial data considered here, the radial reso- 
lution must be at least Nx = 501 to guarantee second- 
order convergence, as the radial features are the domi- 
nant source of numerical error. We note that if we re- 
peat the test above using the same angular grid sizes, 
N(^ = 10, 30, 50, but using instead radial grids with fewer 
points, i.e. N^ = 251,751,1250, the measured conver- 
gence rate is lower, namely n = 1.56. 

For the numerical simulations we present in the re- 
mainder of this article we have chosen grid sizes such 
that the numerical algorithm is always in the second- 
order convergence regime. For more details on the con- 
vergence properties of the radial evolution algorithm, see 

M. 
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FIG. 7: The function x(^) at J^ as a function of Bondi time, 
showing the quasi-normal mode regime oscillations for £ = 1, 
m = 0. The solid line is the output from LEO for A^^ = 11 and 
A''^: = 1001, when the initial data and boundary conditions 
are given as for the convergence test, the dashed line is the 
quasi-normal mode extracted from the data. 



VII. NUMERICAL RESULTS 

A. Quasi-normal Modes 




FIG. 8: Log of the absolute value of the function x(i') s-t J^ 
as a function of Bondi time. Parameters and conditions are 
the same of Fig. [T] The solid line is the output from LEO, 
dashed)line is the quasi-normal mode extracted the data. 

The simulations reported in this and subsequent sec- 
tions are all carried out in compactified, outgoing (re- 
tarded) null coordinates. Because these coordinates al- 
low us to read off scalar radiation patterns at null infinity. 
(The treatment of the inner boundary is as described in 
Sec. I VI CI ) The quality of the waveforms extracted de- 
pends in part on the location of the inner boundary and 
other factors. Wc describe here the method used and 
analyze the sources of error. For the simulations in this 
section we use a grid with sizes Nx = 1500, Nq = 11; 
the initial data corresponds to Eq. ([55]) . with A = 10^""^, 
To = 3M, a = \M, and M = 1. 

To extract the quasi-normal modes we have used the 
free software package Harminv [4^ , which employs a low- 
storage filter diagonalization method (FDM) for finding 
the quasi-normal modes in a given frequency interval. 
This softwarepackage is based on the FDM algorithm 
described in |46l . |47| . The advantage of using Harminv is 
that FDM methods provide better accuracy than what 
can be obtained with a fast Fourier transform (FFT) [0] , 
and are more robust than least-squares fit [ig . We find it 
surprising that this approach, to our knowledge, has not 
been used in the context of reading quasi-normal modes 
in gravitational simulations. 

In performing a fit with Harminv to the scalar field 
waveforms, we find sometimes necessary to factor out, 
at least approximately, the exponential decay of the sig- 
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FIG. 9: The function x at ^ as a function of Bondi time, 
showing the quasi-normal mode regime oscillations for ^ = 2, 
m = 0. The solid line is the output from LEO for A'^^ = 11 and 
Nx = 1001, when the initial data and boundary conditions are 
given as for the convergence test, except that Tin = 2.13M; 
the dashed line is the quasi-normal mode extracted from the 
data. 
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FIG. 10: Log of absolute value of the function x('") at ^ as 
a function of Bondi time. Parameters and conditions are the 
same of Fig. |9] The solid line is the output from LEO, the 
dashed line is the quasi-normal mode extracted from the data. 



nal. This happens when the magnitude of the imagi- 
nary part of the frequency (the decay rate) is compara- 



ble to the real (oscillatory) part, where the FDM method 
fails to find a fitting frequency. In those cases, we pre- 
multiply the signal by an exponentially increasing func- 
tion / = exp(|a;y|f), perform the fit with Harminv, and 
adjust the frequency obtained accordingly. When an an- 
alytic value for the frequency is available, we take its 
imaginary part as the value for ujf. In general, when the 
imaginary part of the frequency is not known, it suffices 
to use a rough estimate of the decay rate, which can ob- 
tained graphically. We also need to decide what range of 
values of u to use to extract this information. We do this 
by plotting the signal x(m) and noting when the waveform 
is clearly periodic with an exponentially decaying enve- 
lope. For example, in Fig.[8l one can clearly see that this 
regime starts at about u = 20M. We take the end of the 
fitting interval when the signal no longer appears to be a 
damped sinusoidal. For initial data of the form ([65|) . with 
i' = 1, we use Harminv to extract the frequency, using as 
the fitting interval u = [20, 70]. The measured frequency 
is w = 0.3076 (5%) - 0.1064 i (9%). Here the values in 
parenthesis indicate the percentage deviation from the 
value calculated in [40] via the WKB method to sixth or- 
der. A comparison of the signal computed and the quasi- 
normal mode fitted is shown in Figs. [7][8l The figures 
show the profiles computed with the three-dimensional 
code (solid line). These profiles are indistinguishable, at 
the resolution of the graph, from the profiles obtained by 
solving numerically the perturbative equation (|44[) for the 
same initial data, thus wc have opted not to show the per- 
turbative solution as is customary. For comparison, we 
have shown instead, in the same graph, the quasi-normal 
mode X = expwu (dashed line) extracted, i.e. the fit pro- 
vided by Harminv. There is some disagreement initially 
between the numerical solution and the fit, as the nu- 
merical solution settles into the dominant quasi-normal 
mode, a process which takes from one to one and half 
cycles of the quasi-normal mode. 



For the same initial data, but with ^ = 2, wc read a 
frequency uj = 0.4971 (3%) - 0.0992 i (2%), in the range 
u = [40, 70], with the comparison between the computed 
signal and [40] shown in Figs. [9lfT0l We have observed 
that the relative percent error for the decay rate is larger 
for £ = 1 because it depends strongly on the value se- 



lected for the location of the boundary. 



Numerical 



experiments with the radial code confirm this and sug- 
gest that, by carefully tuning the location of the inner 
boundary, better accuracy can be achieved for any one 
value of £. Wc have done this only partially in computing 
the frequency for the case £ = 2. Wc want to emphasize 
that the dependence of the frequency on the boundary is 
not a numerical artifact of the code, but a consequence of 
the choice of outgoing null coordinates. This is confirmed 
by numerical experiments with the radial code in ingoing 
coordinates, in which case we find that the frequency can 
be read off with an error of less than 0.1% for the same 
initial data and grid sizes. 
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B. Energy Conservation 
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FIG. 11: Energy conservation as a function of Bondi time for 
l = (S (solid line); 1 = 1 (dotted line); 1 = 2 (long dashed line). 
This calculation was done using the same grid parameters as 
for Fig. [9]except for ri„ = 2.3. For each specific I (line type; 
color) the descending curve corresponds to energy given by 
Eq. (|53fl . The ascending curve corresponds to the algebraic 
sum of Ein = — J Pindu and Eout = j Poutdu. Thus, in 
accordance with Eq. (|56p . the horizontal curve represents the 
global conservation of energy. 

For initial data of the form ^ with ^ = 0, 1, 2, Fig.fTT] 
shows that energy is conserved in the hnear regime. It 
is immediately clear from the graph that the energy con- 
tained on the initial slice is larger the larger the value of 
L In all cases energy is clearly conserved, however, we 
have seen also that if the resolution is not sufhcient for 
a given (., this fact shows up clearly in the graph of en- 
ergy conservation. Thus, we can use energy conservation, 
as well as the results from running the same initial data 
on the radial code, to debug and calibrate the nonlinear 
code, as well as to estimate the evolution time needed 
and its computational requirements. From Fig. II II alone 
the reader might be left to guess as to the extent of the 
deviation of the total energy from a straight line, since 
that deviation is clearly so small that it does not show 
up in the plot for any of three simulations reported in 
Fig. [TTJ Fig. [T^l shows the variation in the energy bal- 
ance AI](w), defined as the percentage variation in S(m) 
relative to the initial value, S](uo), i.e. 



AS = (I](it)/S(0) - 1) X 100, 



(66) 



It can be seen from Fig. [T^ that the relative change 
AE(u) stays below 0.2% during the simulation. We will 
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FIG. 12: Percentage variation in E(?i) with respect to E(0) as 
a function of Bondi time for (. = Q (circles), 1=1 (squares), 
and 1 = 2 (triangles). The graph shows that energy is con- 
served to within less than 0.2% of the energy content of the 
initial surface. 
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FIG. 13: Energy content E{u) as a function of Bondi time 
for: I = Q (circles), 1=1 (squares), 1 = 2 (diamonds), 1 = 3 
(triangles), and £ = 4 (stars). This calculation was done using 
the grid parameters A''^; = 11 and A^^ = 1001. The initial data 
and boundary conditions are the same as in the convergence 
test. 



revisit energy conservation in the context of large resolu- 
tion simulations in Sec. I VII Cl 
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FIG. 14: Energy flow to infinity Eout = J Pout{u)du as a 
function of Bondi time for: £ — (circles), £ = 1 (squares), 
£ = 2 (diamonds), £ — 3 (triangles), and £ — 4 (stars). This 
calculation was done using the same conditions of Fig. 1131 
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FIG. 15: Energy flow into the black hole, Ein = — J Pin{u)du, 
as a function of Bondi time for: £ = (circles), £ = 1 
(squares), £ = 2 (diamonds), £ = 3 (triangles), and £ = 4 
(stars). Both curves for £ = 3 and £ — 4 saturate eventually 
without crossing for u > 30. This calculation was done using 
the same conditions of Fig. 1131 



Fig. [T51 shows the energy content E{u) as a function of 
Bondi time u for a sequence of simulations with initial 
data (p5)) with varying values of £. For lower values of 



^ (^ = 0, 1, 2), the energy content E{u) decays slowly at 
first, then drops rather sharply, and afterwards it decays 
again slowly, at a much lower rate. For higher values of 
£, {£ = 3,4), the energy decays approximately monoton- 
ically from the beginning of the simulation. 

We also observe, see Fig. [T31 that in general, increasing 
£ corresponds to an increase of the energy radiated at J^. 
The oscillations observed in the profiles are higher the 
higher the value of £, as would be expected. The most 
interesting observation in the analysis of energy balance 
arises from Fig. [151 a-nd is the following: for values of £ 
from to 2, the total energy flux towards the black hole 
(as measured by Ein(u) as m — *■ oo) increases with the 
value of £; however, for values of ^ > 2, the total flux of 
energy towards the black hole diminishes with increasing 
values of £. We have confirmed that this is the case with 
the radial code, so this is not a non-linear effect. It is also 
clear that the sudden change of energy for ^ < 2 is due 
to the energy carried away by the scalar field as it falls 
into the black hole. At about i = 2, the radiation into 
the black hole saturates, and for higher values of £, i.e. 
for £ > 2, the centrifugal potential barrier prevents much 
of the field from falling into the black hole. Thus, for 
the same amplitude, configurations with higher angular 
momentum (larger i values) carry more energy, most of 
which will be radiated away and less of which will fall 
into the black hole, so in that sense these configurations 
are proportionally more efficient at carrying energy out 
to J^. 



C. Large resolution simulations 

In order to get a first glimpse of the type of simulations 
that our framework enables us to perform, and to perform 
a final calibration check of the nonlinear code, we select 
initial data given by Eq. ([SS]) . with A = 10~^, tq = 3M, 
a ~ i^M, £ = 8, m = 6, and evolve this configuration un- 
til u = 30Af. This simulation is performed in compacti- 
fied outgoing coordinates, with the treatment of the inner 
boundary as described in Sec. IVICI The plots shown are 
of quantities computed at null infinity, J^. The angu- 
lar grid has size A^^ = 93, that is, there are 372 points 
on a great circle on the sphere, while the radial grid has 
Nj. = 1501 points. This simulation requires 27 hours on 
54 processors, for a total of 1458 processor-hours, or the 
equivalent of two months of a single-processor run. It is 
not by far the largest simulation we could run with our 
framework: we have performed scaling studies that indi- 
cate the code scales linearly well into the 4000+ processor 
range, but it suffices as a demonstration of the resolution 
that can be achieved and the typical turn-around times. 
We assign no particular significance to the initial data 
selected, other than the fact that its angular complex- 
ity provides an excellent test of the code. On any large 
simulation, data analysis and visualization is always a 
challenge. In LEO, visualization is performed by having 
each processor write its own data set at run time, the 
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FIG. 16: Surface plots of x at j?" for u = 2.5M (top) and 
u — 30M (bottom). The parameters of the initial data are 
A = 10"", ro = 3M, a = 0.5M, £ = 8, m = 6. The grid size 
is iVc = 93, N^ = 1501. 



FIG. 17: Surface plots of J J at J" ior u = 2.5M (top) and 
u — 30M (bottom). The parameters of the initial data and 
grid size are the same as in Fig llOl 



individual data files are then post-processed, and graphs 
of the desired quantities generated with Paraview [49|. 
Paraview allows us to easily generate graphs of slices at 
constant coordinate lines and volumetric renderings of 
various fields. Of particular interest to us is the behavior 
of the various metric quantities at null infinity. 

Figs. fT6ll20l display the code variables (j), J, f3, U and W 
as functions on the sphere at null infinity, at u = 2.5M 
and at u = 30 A/. In each case, the graphs show the cor- 
responding field on the six cubed-sphere caps (the gap 
between the caps is the actual size of the spacing be- 
tween grid cells) . The north pole is at the top of the fig- 
ure, rotated 45 degrees towards the observer. For those 
fields that are complex (and which have spin different 
from zero), i.e. J, U, we display for ease of visualization 
the combinations J J and UU, which are real and have 
spin zero. Clearly visible in Fig. \T^ is the m = 6 az- 
imuthal dependence of the field, marked by the presence 
of six maxima and minima. It is also apparent that x 
oscillates in time, as the maxima and minima alternate 
between the top and bottom figures. The angular de- 
pendence is preserved by the evolution, as expected, as 
the only change between the two figures is in the overall 
amplitude (by a factor of w 25 in between the two times 
shown). The graphs of J J, Fig. [T71 are clearly different 



in their angular dependence, showing the presence of var- 
ious harmonics at earlier time. This can be understood 
since in our initial data J{r,x'^) = 0, thus J develops 
from X, i.e. J ~ (S*/*)^, and it is not until later times 
that a definite profile for J has formed. The graph of /?, 
Fig. 1181 shows precisely the angular dependence resulting 
from the contribution f3^r ^ {4>,rY to the source term in 
Eq. (|32p . and remains constant throughout the simula- 
tion, up to an overall amplitude. The graphs of U and 
W , Figs. [TM^ show higher order angular dependence 
arising from the angular derivatives of C/, which in turn 
are essentially driven by the source term in the Eq. (|34p . 
i.e. hy U ^ Q ^ <f)d<j). 

Fig. [5T] shows again that energy is conserved during 
the entire simulation. The variation in the energy balance 
S(w) is well below 1%, thus S(u) is indistinguishable from 
a straight line at the resolution of the graph, as noted 
in Sec. IVII Bl To more fully appreciate to what extent 
energy is conserved, the graph insert in Fig. [21] shows 
the percentage variation of S(u), normalized to its value 
at M = 0. From the graph insert we see that the total 
energy varies by at most 0.025% during the simulation. 
As we stated earlier, energy conservation is a requisite for 
accuracy in the waveforms; note that energy is conserved 
during this simulation to within tighter limits than in the 
simulations of Sec. IVII Al which is due primarily to the 
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FIG. 18: Surface plots oi f3 a.t J" for u = 2.5M (top) and 
u = 30M (bottom). The parameters of the initial data and 
grid size are the same as in Fig |16l 



FIG. 19: Surface plots of UU at ^ for u = 2.5M (top) and 
u — 30M (bottom). The parameters of the initial data and 
grid size are the same as in Fig |16l 



increased angular resolution. For this simulation we took 
N(^ ~ 93 in order to accurately resolve the higher order 
harmonic angular dependence, where in Sec. IVII Al we 
set N(^ = 11; this amounts to (approximately) an 8-fold 
increase in resolution, and for the same radial resolution, 
a 64-fold increase on the computing resources required. 

We would be remiss if we did not discuss at least briefly 
the performance characteristics of our code. As part of 
our calibration and testing, we have performed detailed 
profiling studies, which we will not go into detail here. 
Suffice to say that in its current configuration, the code 
performs at approximately 20% of peak on the Cray XT3. 
Its weak scaling is linear (that is, its performance solv- 
ing progressively larger configurations, while keeping the 
load per processor constant, scales linearly with the num- 
ber of cores), while running on up to 4056 CPU cores on 
the Cray XT3 at PSC. 



VIII. CONCLUDING REMARKS AND 
OUTLOOK OF FUTURE WORK 

We have presented a new computational framework 
(LEO), which we can use to perform large-scale, high- 
resolution calculations in the context of the charac- 
teristic approach in numerical relativity. This highly 



parallel and easily extensible implementation has been 
used to solve the model problem of a massless scalar 
field minimally coupled to gravity (the three-dimensional 
Einstein-Klein-Gordon problem). We have shown that 
the nonlinear code thus implemented is globally second- 
order convergent, and how accurately we can follow 
quasi-normal mode ringing. We have studied the balance 
of energy for a number of initial data sets with different 
angular structure. Aside from the interesting result of 
energy flow saturation through the Schwarzschild hori- 
zon, the LEO framework offers a good prospect to study 
new configurations beyond the linear regime and the grid 
sizes used in this work. 

Future directions we are currently exploring include 
the application of the LEO framework to a consistent, 
quasilinear, fully first-order formalism derived from p^ . 
and the extension of the model problem considered here 
to massive scalar fields. The later case is particularly 
important because it will allow us to simulate a boson 
star orbiting a black hole. We will compare the per- 
formance of the characteristic framework in ingoing ver- 
sus outgoing null coordinates in the extraction of quasi- 
normal modes, and in the study of nonlinear effects in the 
neighborhood of the central black hole. The underlying 
framework can be applied equally well to 3 -I- 1 formula- 
tions of the Einstein equations in spherical coordinates. 
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FIG. 20: Surface plots of W at J^ for u = 2.5M (top) and 
u — 30M (bottom). The parameters of the initial data and 
grid size are the same as in Fig |16l 



stantial improvements in the computed waveforms, and 
this issue remains to be addressed. Although of potential 
importance for the accuracy of gravitational waveforms, 
such a study lies outside the scope of the present work. 

2.5 




30 



FIG. 21; Energy conservation for the simulation which gen- 
erated the results shown in Figs. I16H2UI The solid line cor- 
responds to the energy content E{u) at successive times, the 
dashed line to the sum of the energy radiated through the 
inner (Ein{u)) and outer {Eout{u)) boundaries, and the dot- 
dashed line to the sum, E(u) = E{u) + Ei„{u) + Eout{u), 
respectively. The insert graph shows the percentage variation 
in E(u), relative to its final value at u = 30. 



in particular to a generalization to three dimensions of 
the Bondi-Sachs gauge of [4^, and finall y t o matched 
3-1-1 and characteristic evolutions [42, [50, l5l| . 

We have not addressed in the present work some out- 
standing problems with the calculation of the News [301 1 
some of which arise from second angular derivatives of 
the metric fields at ^ entering in the computation of the 
News, a feature which can lead to substantial propaga- 
tion of errors. We have observed during our simulations 
that the metric fields computed at J^ are smooth, as 
evident in Figs. [16ll20l (and so are those fields which rep- 
resent their angular derivatives, although these are not 
shown here). It should be noted that, in our formula- 
tion [^ , the first angular derivatives of some fields have 
been promoted to auxiliary variables, for which a hyper- 
surface equation is integrated radially. (We do this for 
those fields whose second angular derivatives enter in the 
computation of the News). In practice, this means that 
only first order angular derivatives of any of the fields 
we evolve need to be computed to calculate the news. It 
is possible that the cubed-sphere approach leads to sub- 
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